#############################################################
#########    Altenrtive Specification ############
##################   Cubic Polynomial        #############
###############      October 10, 2018    ##################
#######  Check rerun: December 16, 2022 #######

rm(list=ls())
library(foreign)
library(plyr)
library(multiwayvcov)
library(sandwich)
library(lmtest)
library(stargazer)
library(ggplot2)


data=read.csv("~/Dropbox/Personal Research 2017/replications/karn_nov16.csv")
names(data)
summary(data$Latitude)
summary(data$Longitude)
summary(data$border1)
summary(data$border2)
summary(data$Slope)


#######################
##### Distances #####

#Distance to Mysore-Bombay Border
rd10.mb=data[which(data$NEAR_DIST_border1<10000),] #20 km


table(rd10.mb$border1)

#Distance to Hyderabad-Bombay Border
rd10.hb=data[which(data$NEAR_DIST_border2<10000),] #20 km

table(rd10.hb$border2)




# mysore 
rd10.mb$cubic=polym(rd10.mb$Latitude, rd10.mb$Longitude, degree=3, raw=TRUE) #cubic poly
summary(rd10.mb$cubic)

health.mys=lm(health_binary~border1+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+rd10.mb$cubic, data=rd10.mb)
summary(health.mys)
health.mys.cl=cluster.vcov(health.mys, rd10.mb$dist_name)
health.mys.se=sqrt(diag(health.mys.cl)) #cluster standard errors


pucca.mys=lm(pucca_binary~border1+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+rd10.mb$cubic, data=rd10.mb)
summary(pucca.mys)
pucca.mys.cl=cluster.vcov(pucca.mys, rd10.mb$dist_name)
pucca.mys.se=sqrt(diag(pucca.mys.cl)) #cluster standard errors


# 
# stargazer(mcl2, mcl3, digits=3,
#           column.labels = c("Health Facilities", 
#                             "Pucca Roads"))
# 


# hyderabad 
rd10.hb$cubic=polym(rd10.hb$Latitude, rd10.hb$Longitude, degree=3, raw=TRUE) #cubic poly
summary(rd10.hb$cubic)

health.hyd=lm(health_binary~border2+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+rd10.hb$cubic, data=rd10.hb)
summary(health.hyd)
health.hyd.cl=cluster.vcov(health.hyd, rd10.hb$dist_name)
health.hyd.se=sqrt(diag(health.hyd.cl)) #cluster standard errors

pucca.hyd=lm(pucca_binary~border2+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+rd10.hb$cubic, data=rd10.hb)
summary(pucca.hyd)
pucca.hyd.cl=cluster.vcov(pucca.hyd, rd10.hb$dist_name)
pucca.hyd.se=sqrt(diag(pucca.hyd.cl)) #cluster standard errors



# stargazer(mcl2, mcl3, digits=3,
#           column.labels = c("Health Facilities", 
#                             "Pucca Roads"))




stargazer(health.mys, pucca.mys, health.hyd, pucca.hyd, se=list(health.mys.se, pucca.mys.se, health.hyd.se, pucca.hyd.se), digits=3, 
          omit=c("TOT_POP", "TOT_SC", "TOT_ST", "Slope", "TerrainRug", 
                 "cubic1.0", "cubic2.0", "cubic0.1", "cubic1.1",
                 "cubic2.1", "cubic0.2", "cubic1.2", "cubic0.3", "cubic3.0"), 
          dep.var.labels=c("Health Centers", "Paved Roads", "Health Centers", "Paved Roads"), 
          covariate.labels = c("Indirect Rule (Mysore)", "Indirect Rule (Hyderabad)", "Constant"),
          add.lines = list(c("Controls","\\checkmark","\\checkmark","\\checkmark","\\checkmark")), 
          omit.stat = c("rsq", "f", "adj.rsq", "ser"))


